*     logical function returns .TRUE. if CZ lies within the 4 nodes in
*     the array CZNODE
*     if CZ is closer to any of the sides (with half lengths stored in
*     the array RHL) than the tolerance RTOL the CZN returns new
*     coordinates for CZ that are RTOL away from that side
      LOGICAL FUNCTION LFSAIN(CZ, CZNODE, RHL, RTOL, CZN)
      IMPLICIT NONE

      COMPLEX*16 CZ, CZNODE(*), CZN
      DOUBLE PRECISION RHL(*), RTOL

      INTEGER IN, I2, I3
      COMPLEX*16 CZ1, CZ2, CZL, CZL2, CZ0, CZLN, CZ3
      DOUBLE PRECISION RLHL, RLHL2, RX, RY, RY2, RTH
      COMPLEX*16 CDLOG

      LFSAIN=.FALSE.
      RTH=0.
      CZN=CZ
*     loop over all sides to check if CZ is RTOL away from sides and if
*     not then calculate new coordinates CZN
*     RETURN statement after action in case CZ at a corner
      DO IN=1,4
*        the element is translated and rotated such that the center is at
*        the origin and the element along the x-axis with endpoint 1 at x1=
*        -.5*(length of element), endpoint 2 at x2=.5*(length of element)
         I2=IN+1
         IF (IN.EQ.4) I2=1
         CZ1=CZNODE(IN)
         CZ2=CZNODE(I2)
         RLHL=RHL(IN)
         CZL=(2.*CZN-(CZ2+CZ1))*RLHL / (CZ2-CZ1)
         RX=REAL(CZL)
         RY=AIMAG(CZL)
         IF (ABS(RY).LT.RTOL) THEN
*        check left corner for first element
            IF (IN.EQ.1) THEN
*              check with previous element for left corner
               CZ0=CZNODE(4)
               RLHL2=RHL(4)
               CZL2=(2.*CZN-(CZ1+CZ0))*RLHL2 / (CZ1-CZ0)
               RY2=AIMAG(CZL2)
               IF (ABS(RY2).LT.RTOL) THEN
*                 indeed at corner
                  CZL2=(2.*CZ0-(CZ2+CZ1))*RLHL / (CZ2-CZ1)
*                 only action at corner if the two elements are
*                 not (almost) parallel
                  IF (ABS(ATAN2(AIMAG(CZL2), REAL(CZL2)+RLHL))
     &               .GT.RTOL) THEN
                     IF (RY.GE.0..AND.RY2.GE.0.) LFSAIN=.TRUE.
                     IF ((AIMAG(CZL2).LT.0).AND.
     &                  (RY.GE.0..OR.RY2.GE.0.)) LFSAIN=.TRUE.
                     IF (LFSAIN) THEN
*                       inside
                        CZLN=
     &              -RLHL+RTOL/AIMAG(CZL2)*(CZL2+RLHL+2.*RLHL2)
                     ELSE
*                       outside
                        CZLN=
     &              -RLHL-RTOL/AIMAG(CZL2)*(CZL2+RLHL+2.*RLHL2)
                     END IF
                     CZN=.5*(CZLN*(CZ2-CZ1)/RLHL+CZ2+CZ1)
                     RETURN
                  END IF
               END IF
            END IF
*           check right corner if not checking the last element
            IF (IN.NE.4) THEN
*              check with next element
               I3=I2+1
               IF (I2.EQ.4) I3=1
               CZ3=CZNODE(I3)
               RLHL2=RHL(I2)
               CZL2=(2.*CZN-(CZ3+CZ2))*RLHL2 / (CZ3-CZ2)
               RY2=AIMAG(CZL2)
               IF (ABS(RY2).LT.RTOL) THEN
*                 indeed at corner
                  CZL2=(2.*CZ3-(CZ2+CZ1))*RLHL / (CZ2-CZ1)
*                 only action at corner if the two elements are not
*                 (almost) parallel
                  IF (ABS(ATAN2(AIMAG(CZL2), REAL(CZL2)-RLHL))
     &          .GT.RTOL) THEN
                  IF (RY.GE.0..AND.RY2.GE.0.) LFSAIN=.TRUE.
                  IF ((AIMAG(CZL2).LT.0.).AND.
     &            (RY.GE.0..OR.RY2.GE.0.)) LFSAIN=.TRUE.
                  IF (LFSAIN) THEN
*                    inside
                     CZLN=
     &              RLHL+RTOL/AIMAG(CZL2)*(CZL2-RLHL-2.*RLHL2)
                  ELSE
*                    outside
                     CZLN=RLHL-RTOL/AIMAG(CZL2)*(CZL2-RLHL-2.*RLHL2)
                  END IF
                  CZN=.5*(CZLN*(CZ2-CZ1)/RLHL+CZ2+CZ1)
                  RETURN
               END IF
            END IF
          END IF
*     check if along element
          IF (ABS(RX).LT.RLHL) THEN
            IF (RY.LT.0.) THEN
              RY=-RTOL
            ELSE
              RY=RTOL
            END IF
            CZN=.5*(CMPLX(RX, RY)*(CZ2-CZ1)/RLHL+CZ2+CZ1)
          END IF
        END IF
        RTH=RTH+AIMAG(CDLOG((CZN-CZ2)/(CZN-CZ1)))
      END DO
      IF (RTH.GT. 3.14) LFSAIN=.TRUE.
      RETURN
      END
